function [Rs, Rs2, cm, tau1] = Rsestimate()
%
% Rsestimate computes the series resistance, cell cap, and a time contstant (fastest) for 
% a voltage clamp step response. Based on equations in Joe Santos-Sacchi's 1993 paper
% (Voltage-dependent ionic conductances of type I spiral ganglion cells from the guinea pig inner ear.
% J Neurosci. 1993 Aug;13(8):3599-611. )
% Initial version, 9/1/06 P. Manis

% Modified 9/8/06 Paul B. Manis, Ph.D.
%
% Now using FIT parameters to estimate Q (integral of fitted exponential),
% and the 3 currents (A1, A2 associated with a double exponential fit) and
% A0, the "steady-state" due to Rm alone.
%
% Tests show that the estimates are more accurate this way, and Rs2 is now
% returned, corresponding to the fast component only.
% Useful Reference: Jackson, M.B., Cable analysis with whole-cell patch
% clamp, Biophys. Journal. 61: 756, 1992.

% Errors in estimates of Rs can still arise from:
% 1. Sampling too slowly
% 2. duration of voltage step too short to get good estimate of second
% exponential 
% 3. Noise
% 4. extreme parameters (Rs > 0.5*Rm).
% 5. Incorrect compensation for the pipette capacitance (transmural Cp)
%



% define parameters of the test pulse
vc = -10; % voltage step, in mV
tstart = 45.6; % start time of the test pulse, in msec
tend = 75.6; % end time of test pulse in msec
testmode = 0; % flag to control test mode versus real data mode
%
% in test mode, we evaluate the results of the algorithm over a range of input parameters.
if(testmode)
    tstep = 0.05; % time step (in test mode only; real data mode uses sample time step from data array)
    offcurr = -300;
    x=0:tstep:150;

    rst = [0.005 0.010 0.020 0.050 0.1]; % Gohm
    rint = [2 1 0.5 0.2 0.1]; % Gohm
    cmt = [5 10 50 100 200]; % pF
%     rst = [ 0.010 0.050]; % Gohm
%     rint = [1 0.5  0.1]; % Gohm
%     cmt = [5 20 100]; % pF
    secondtau = 0.1;
    xst=find(x>=tstart);
    xend=find(x>=tend);
    xst = xst(1);
    xend = xend(1);
    for ir = 1:length(rst)
        for ii = 1:length(rint)
            for ic = 1:length(cmt)
                    y=zeros(length(x),1) + offcurr;
                    y(xst:xend) =  y(xst:xend) + (vc/rint(ii)) + (vc/rst(ir))*exp(-(x(xst:xend)-tstart)/(rst(ir)*cmt(ic)))';
                    y(xst:xend) = y(xst:xend) + (secondtau*vc/rst(ir))*exp(-(x(xst:xend)-tstart)/15)';
                    [Rs, Rs2, Rin cm, tau1, a2, tau2] = santosrs(x, y', vc, tstart, tend, tstep);
                fprintf(1, 'Rs=%8.4f Gohm (%8.4f) [Rs2 = %8.4f]  Rin = %8.4f Gohm (%8.4f)  Cm=%7.2f pF (%7.2f)  tau1 = %8.3f msec (%8.3f)  [a2 = %8.3f  tau2 = %8.3f]\n', ...
                    Rs, rst(ir), Rs2, Rin, rint(ii), cm, cmt(ic), tau1, ((rst(ir)*rint(ii))/(rst(ir)+rint(ii)))*cmt(ic), a2, tau2 );
            end;
        end;
    end;
else
    % evaluate real data here
    
    a=load('/Users/pmanis/Desktop/06414004.txt', 'ascii');
    % get and scale the time array
    x=a(:,1); % time base
    x=x*1000; % scale for sample rate
    tstep = x(2)-x(1);
    y=mean(a(:,2:end)');  % get mean of the current trace

    Rs = []; Rs2 = [];
    cm=[]; tau1=[]; Rin=[];
    for i = 2:size(a,2)
        y = a(:,i);
        [Rs(i), Rs2(i), Rin(i), cm(i), tau1(i), a2(i), tau2(i)] = santosrs(x, y, vc, tstart, tend, tstep);
        fprintf(1, 'Rs=%8.4f Gohm  [Rs2 = %8.4f] Rin = %8.4f Gohm   Cm=%7.2f pF   tau1 = %8.3f msec a2 = %8.3f  tau2 = %8.3f\n', ...
            Rs(i), Rs2(i), Rin(i),  cm(i), tau1(i), a2(i), tau2(i) );
    end;
    h = findobj('tag', 'Rssummary');
    if(isempty(h))
        h = figure('tag', 'Rssummary');
    else
        figure(h);
        clf;
    end;
    xi = 1:length(Rs);
    subplot('Position', [0.1, 0.55, 0.8, 0.35]);
    plot(xi, Rs2*1000, 'rx');
    ylabel('Rs (Mohm)');
    subplot('Position', [0.1, 0.30, 0.8, 0.15]);
    plot(xi, cm, 'g+');
    ylabel('Cm (pF)');

    subplot('Position', [0.1, 0.10, 0.8, 0.15]);
    plot(xi, tau1, 'bs');
    ylabel('Tau (ms)');


end;


function [Rs, Rs2, Rin cm, tau1, a2, tau2] = santosrs(x, y, vc, tstart, tend, tstep)

% function to do the real work
xst=find(x>=tstart);
xend=find(x>=tend);
xst = xst(1);
xend = xend(1);

h = findobj('tag', 'Rstrace');
if(isempty(h))
    h = figure('tag', 'Rstrace');
else
    figure(h);
    clf;
end;



xlim([0, 30]);
i0 = mean(y(1:xst-1));
i1 = mean(y((xend-5):xend));

Rin = vc/(i1-i0); % and Rin = Rm + Rs... 
dt=x(2)-x(1);

% get the indices corresponding to the step and response we need to measure
% and find out where the peak value is located


[Ipeak, jpeak] = min(y(xst:xend));
xfitstart = xst+jpeak-1;
xfit = x(xfitstart:xend)-tstart;
yfitd = y(xfitstart:xend)-i0; % to simplify life, subtract the main offset
plot(xfit, yfitd, 'ro');
[Ipeak, jpeak] = min(yfitd); % get new values...
i1 = mean(yfitd(floor(end*0.8):end));
Imax = abs((Ipeak-i1)*2);

% calculate total charge: summation is for a trapezoidal integration method
Q = dt*abs(sum(yfitd(2:end-1)-i1)+0.5*(yfitd(1)-i1)+0.5*(yfitd(end)-i1));

% generate estimates for the exponential curve fit

a0 = i1;
a1 = Ipeak-i1; % put all the amplitude in the fast component first
t1a = 0.5;

[im, jm] = min(((yfitd-i1)-0.2*min(yfitd-i1)).^2); % find where it drops to 20% of the peak value
if(jm <= 2)
    jm = 5; % fit at least 5 points
end;

% perform a single exp fit on the initial part of the curve first to optimize parameter estimates
% Note: we allow the amplitude and tau to vary, but not the baseline
% constant (VP = [0 1 1])
[FP1, CHISQ1, NITER1] = mrqfit('exponential', ...
     [a0, a1, t1a], xfit(1:jm), yfitd(1:jm), ...
     [], [0, 1, 1], [-5000, -Imax, tstep/2], [5000, Imax, 25], 400, 1e-6 );

%
% use the single exp fit to seed the fast component, then fit again with a
% double exp fit.

a2 = a1*0.1; % scale to fast amplitude
t2a = 10;

[FP2, CHISQ2, NITER2] = mrqfit('exponential', ...
    [FP1(1), FP1(2), FP1(3), a2, t2a ], xfit, yfitd, ...
    [], [1 1 1 1 1], [-5000, -5000, FP1(3)/2, -Imax, FP1(3)*3], [5000, 0, FP1(3)*2, 0, 50], 500, 1e-7 );

tau1 = FP2(3);
tau2 = FP2(5);
a2 = FP2(4);
hold on
    yfit2 = FP2(1) + FP2(2)*exp(-xfit/tau1) + FP2(4)*exp(-xfit/tau2);
    sserr2 = sqrt(sum((yfit2-yfitd).^2))/length(yfitd);
    yfit1 = FP1(1)+ FP1(2)*exp(-xfit/FP1(3));
    sserr1 = sqrt(sum((yfit1-yfitd).^2))/length(yfitd);
%    fprintf(1, 'Fit info: sserr1: %f   sserr2:  %f   ratio: %f\n', ...
%        sserr1,  sserr2, sserr2/sserr1);
if(sserr2/sserr1 < 0.8)
    % now, if we have 2 exps, we need to revise Rin and Q so that Rin
    % reflects both the steady state and the slow (dendritic) charging
    % and Q reflects ONLY the charge associated with the fast (somatic)
    % component.
    Rin = vc/(FP2(1)); % Rs is only due to the fast component.
    Q = abs(FP2(2)*FP2(3)); % int [0->inf] of Ae(-bt) dt = A/b; of Ae(-t/tau) = A*tau.
%    fprintf(1, '2exp: Rs2 = %8.4f   ', vc/FP2(2));
plot(xfit, yfit2, 'g-');
else
    Rin = vc/(FP2(1)); % Rs is only due to the fast component.
    Q = abs(FP2(2)*FP2(3)); % get from fit
    plot(xfit, yfit1, 'b-');
    tau1 = FP1(3);
end;

vca = abs(vc);
Rs2 = vc/FP1(2);

% fprintf(1, '1exp: Rs2 = %8.4f   ', vc/FP1(2));

% calculate the resulting values using equations from JSS paper
% verified that the equations are algebraically correct (P. Manis, 9/5/06).
Rs = (Rin*tau1*vca)/(Q*Rin+tau1*vca);
Rm = Rin - Rs;

cm = (Rin^2*Q)/((Rm^2)*vca);

return;
